Finite Temperature Dynamical Structure Factor of the Heisenberg— Ising Chain 
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We consider the spin-1/2 Heisenberg XXZ chain in the regime of large Ising-hke anisotropy A. 
By a combination of duahty and Jordan-Wigner transformations we derive a mapping to weakly 
interacting spinless fermions, which represent domain walls between the two degenerate ground 
states. We develop a perturbative expansion in 1/A for the transverse dynamical spin structure 
factor at finite temperatures and in an applied transverse magnetic field. We present a unified 
description for both the low-energy temperature-activated response and the temperature evolution 
of the T=0 two-spinon continuum. We find that the two-spinon continuum narrows in energy 
with increasing temperature. At the same time spectral weight is transferred from the two-spinon 
continuum to the low energy intraband scattering continuum, which is strongly peaked around the 
position of the (single) spinon dispersion ("Villain mode"). 

PACS numbers: 75.10.Jm, 75.10.Pq, 75.40.Gb 



I. INTRODUCTION 

The spin-1/2 Heisenberg XXZ chain is a paradigm in 
low-dimensional quantum magnetism. Its Hamiltonian 



+1" 



qy qy 



(1) 



where h is an external magnetic field and A controls 
the exchange anisotropy. The one-dimensional case is 
particularly significant because for a field h parallel to 
the z direction, the Hamiltonian is integrable and the 
spectrum of the spin chain may be extracted exactly. 
Furthermore the Hamiltonian (1) is thought to provide 
a realistic description of several quasi-lD experimen- 
tal compounds. Examples include CS2C0CI4 for which 
A = 0.25,^^ CsCoBrs,^'^ CsCoCla^^ and TlCoCla^^ all 
of which have A ~ 7. 

For A > 1, = and T = the XXZ spin chain is 
in a Neel phase. The fundamental excitations take the 
form of gapped fractionalized spin-1/2 quantum solitons 
known as spinous. Strictly in the Ising limit, A = 00, 
the spinous can be identified simply as domain walls 
with gap A/2. Experiments have established the exis- 
tence of several manifestations of the XXZ model in the 
Ising regime. In some cases these experiments have 
also probed the effects of temperature^' and transverse 
field^^ on dynamical spin-spin correlations. 

The measure of dynamical correlations known as the 
dynamical structure factor is an important quantity in 
the study of quantum magnets. This is for two reasons: 
firstly, the dynamical structure factor is directly mea- 
surable by inelastic neutron scattering experiments and 
secondly the nature of the dynamical response is highly 
specific to the system in question, so that it serves as 
a characterization tool. A particular feature that one 
would like to understand in the case of the XXZ chain is 
the finite temperature low energy spin response known as 



the 'Villain mode'.^^ This response due to scattering be- 
tween domain wall pair states has been observed in Refs. 
[9,13]. As this response occurs only at finite tempera- 
tures it necessitates a theory that accounts for thermal 
fluctuations. A recent analysis of the continuum limit of 
two gapped integrable quantum spin chains^^ has shown 
that at raised temperatures the effect of thermal fluctu- 
ations cannot be described in terms of a simple thermal 
decoherence or relaxation time picture. Instead markedly 
asymmetric thermal broadening of single particle modes 
is observed. This paradigm has been found to be in agree- 
ment with theoretical^^' and experimental^^ studies of 
the spin-1/2 Heisenberg chain with strongly alternating 
exchange, a model which is gapped but not integrable. 
In contrast the gapped excitations in the spin-1/2 XXZ 
chain occur only in pairs; it is then interesting to try and 
understand the thermal evolution of the resulting two 
particle response in addition to that of the Villain mode. 

Despite the advantages afforded by integrability, the 
task of calculating correlation functions (and hence dy- 
namics) for the XXZ chain is still far from simple. 
First order perturbative treatments around the A = 00 
limit^'^^'^^ and 1/S expansions'^ have given some insight. 
In recent years significant progress has been made for the 
A > 1 regime, both via Bethe's Ansatz'^ and a different 
exact technique which works directly with the thermo- 
dynamic limit.'' This has lead to an exact expression 
for the transverse dynamical correlations at T = 0.21.23 
Results at finite temperature are generally still limited 
to asymptotically exact thermodynamic quantities'^ 
However there have been promising advances that rely on 
generalizing multiple integral representations for time de- 
pendent correlation functions to finite temperatures.'^^' 
Currently these methods have not yielded expressions for 
the most experimentally relevant quantity, the dynamical 
structure factor. 

In this paper we present a perturbative calculation 
for the transverse spin response in the A ^ 1 limit, 
valid at finite temperatures and capable of incorporat- 
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ing a transverse field. We note that the temperature 
dependence of the dynamical structure factor in the crit- 
ical — 1 < A < 1 regime has been determined by ex- 
act diagonalization of short chains and very recently by 
Quantum Monte Carlo and DMRG computations.^^ The 
usual perturbative approach to the XXZ chain uses the 
Jordan- Wigner transformation^^ to produce an expan- 
sion in powers of A. This is suitable for investigating 
the XY, |A| <C 1 case but inadequate here. Instead we 
first perform a Kramers- Wannier^^ duality transforma- 
tion to a new Hamiltonian in terms of dual operators. 
A Jordan- Wigner transformation from these dual oper- 
ators to spinless fermions then leads to a controlled ex- 
pansion in 1/A. An equivalent mapping has been used 
previously, coupled with mean field theory, to find the ap- 
proximate excitation spectrum in the Ising phase. We 
take an alternative approach resumming certain terms in 
the expansion to all orders, in order to take account of 
both quantum and thermal fluctuations. 

The structure of the paper is as follows. First in Sec. 
II we discuss symmetries of the Hamiltonian and the dy- 
namical structure factor. Second, in Sec. Ill we trans- 
form the Hamiltonian into a form suitable for the expan- 
sion. In Sec. IV we describe the perturbative expansion 
of the transverse spin-spin correlator. In Sec. V we ex- 
plain how to resum certain terms in this expansion in 
order to obtain finite results. In Sec. VI we discuss the 
behaviour of the dynamical structure factor for a range 
of parameters. Sec. VIII contains some brief concluding 
remarks. 



II. SYMMETRIES OF THE HAMILTONIAN 
AND THE STRUCTURE FACTOR 

We now describe the symmetries of the Hamiltonian 
and their consequences for spin-spin correlation func- 
tions. For h = the Hamiltonian, Eq. (1), is invariant 
under arbitrary rotations 7l^{(t)) around the z-axis as well 
as under rotations by tt around the x axis, 7^^(7r), which 
entail the mapping 

~^ ' 

qy,z qy,z 
^ ' 

The two types of symmetry operations do not commute, 
but we can diagonalize the Hamiltonian simultaneously 
with either or with the generator IZ^ij:) of the Z2 
symmetry. This in turn implies that all off-diagonal 
spin correlators vanish for T > 0, by the following argu- 
ments. When considering the thermal expectation values 
{S^S^) with a = x^y we choose a basis of energy eigen- 
states in which the total is diagonal. Then carrying 
out a rotation by tt around the z-axis sends —^n^ 
a = x^y and as a result 

{S:S^)=-{S^S^J = 0, a = x,y. (2) 

On the other hand, when considering (6'^ 6'^) we use a 
basis of simultaneous eigenstates of H and IZ^ (tt) to carry 



out the thermal trace. Under the Z2 symmetry the ther- 
mal expectation value is negated, leading to 

{s;isyj = -{s:syj = o. (3) 

This shows that in the absence of the transverse magnetic 
field all off-diagonal elements of the dynamical structure 
factor vanish. In the presence of a finite transverse field 
we only have the Z2 to work with. Concomitantly for 
/i > one finds {S^Syj = {S^S^) = but is no 

longer required to vanish by symmetry and as a result 
acquires a finite value. 



III. TRANSFORMATIONS OF THE 
HAMILTONIAN 

In order to proceed we aim to re-express the Hamilto- 
nian in terms of spinless fermions by means of a Jordan- 
Wigner transformation in such a way that we analyze 
the resulting interacting fermion Hamiltonian by stan- 
dard perturbative methods. In order to achieve this, 
the AJ^^ 'Sn'^n-\-i P^^^ of the Hamiltonian (1) must be 
mapped to an expression quadratic in fermions. One ap- 
proach for doing this is outlined in Appendix A, another 
one is discussed in detail next. In the following we con- 
sider a transverse magnetic field applied along the x di- 
rection. We note that the effects of of a transverse field 
in the critical region of the XXZ-chain — 1 < A < 1 have 
been studied in some detail in Refs 37. 



A. Duality Transformation 

We work on the infinite chain so that we can ignore 
boundary effects. The Hamiltonian, Eq. (1), in terms of 
Pauli spin matrices is given by 

H — Ha + Hh, 

= E «+i + 1 («+i + ^ 

n n 

Hu = \Y.<- (4) 

n 

The Kramers- Wannier duality transformation is defined 

by 

J<n+1 

^n+1/2 ^ ~^Mn+l/2Mn+l/2- (^) 

This transformation defines operators on a dual lattice 
and maps an Ising chain from its ordered to its disor- 
dered phase. Applying the transformation to the XXZ 
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Hamiltonian we find 



with raising and lowering operators 



+ 1/2 



+ /^n-l/2Mn+3/2 Mn+l/2Mn-l/2Mn+3/2) ' 
= y^Mn-l/2Mn+l/2- (6) 



(9) 



and then use the Jordan-Wigner transformation: 



We note that applying the duality transformation to a 
finite open chain leads to a dual Hamiltonian containing 
additional boundary terms that in particular ensure a 
doubly degenerate ground state in the thermodynamic 
limit. This is most easily seen for A ^ h ^ 1, i.e. 
the Ising model in a transverse field. In this limit the 
mapping gives 



N-l 



N 



n—l n—1 
N-l 

-^^Yl /^n+l/2Mn+3/2 + ^ ^n+1/2 + ^3/2 ' (7) 



N-l 



n=l 



Now as h/A the two ground states of the origi- 
nal Hamiltonian become the familiar Neel states, (cr^) = 
~{^n-\-i)' contrast, for the dual Hamiltonian the 
ground state is given by (^^+1/2) — ~^ ^^^^ n < N. 
The twofold degeneracy is then maintained by the TVth 
dual spin, (jln +1/2 which is free to point in either direc- 
tion. 

In the following we will be interested only in bulk cor- 
relations of operators that are local under the duality 
transformation. Hence the boundary terms do not play 
a role and will be dropped. 



(10) 



1. Spin Operators 



We first consider the transformation of the lattice spin 
operators under the mappings. Crucially, the transverse 
spin operator is local under the transformations 



On the other hand. 



both 



t t 

^n—l^n 

and 



h.c. 



(11) 



(7^ acquire Jordan- 
Wigner strings. As a result, our formalism will allow us to 
determine the xx-component of the dynamical structure 
factor only. It follows from the symmetry considerations 
above that in absence of a transverse magnetic field this 
suffices to determine all transverse correlations. 



B. Fermionic Representation 



In order to proceed further it is necessary to map the 
spins to fermions. We perform a rotation of spin axes 



^n+1/2 ~^ ^n^ ^n+1/2 



(8) 



2. Hamiltonian 



After the Jordan-Wigner transformation the Hamilto- 
nian takes the form 



i^A = ^ A ^ cicn + ^ [4-lCn+l - " 4-l4cnCn+l + ^-l^Cncl+l + ^.c] + COUSt. (12) 



and 

^h = \Y. - ^n-l^n + ^-C-) • (13) 

n 

We now write the Hamiltonian as a sum of two pieces, 
H = H2-\- containing terms quadratic and quartic in 



the fermionic operators respectively. The quartic (inter- 
action) terms are O(A^) while the quadratic pieces mix 
orders 0{A) and O(A^). The external field only appears 
in the quadratic part of the Hamiltonian. After taking 
the Fourier transform of the quadratic part, i^2, of the 
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Hamiltonian we find 



^2 = ^E(4 c-,) (. 



Ak iBk 
-iBk -Akj \c 



(14) 



where we have defined the new functions 



with Ak = 2A+4 cos(2A;) + ^ cos(A:) and Bk = ^ sin(2A:) + 
^sin(/c). This can be diagonalized by a Bogoliubov 
transformation of the form 



icos{Ok) -sm{Ok)\ f OLk 
sin(6>/c) -icos{Ok)) \a^_. 



so that 



(15) 



(16) 



f{kiM){k3M) = cos(/ci - k^) - cos(/c2 - k^) 

- cos(/ci - /cs) + cos(/c2 - /ca), 

dikiMM) = sin(A;i - A;2) + sin(A;2 - ks) 
+ sin(A:3 - /ci). 



(19) 



(20) 



with tan(2(9/e) = ^/c/A/, and Ck = ^a/A^ +5^. 

Now we consider the quartic part of the Hamiltonian: 

n 

(17) 

Taking the Fourier transform and manipulating indices 
leads to 



Ha — - ^2 ^fcl+fc2+/C3+/C4 

fcl ,/C2 ,^3,^4 

X [/(fcife2)(fe3/C4)4i42^-fc3C-/C4 

2 

+ 3 



(^^(/ci/C2/C3)4i42 43^-^4 +h.C.) ], (18) 



The function / is antisymmetric under exchange of the 
two momenta within a pair of brackets. The function 
g is symmetric for cyclic permutations of its momen- 
tum arguments and antisymmetric otherwise. For exam- 
ple /(fei,fe2)(/C3,/C4) = ~/(/C2,fei)(fe3,fc4) = f {ka M) {ki M) 

dikiMM) = 9{k2MM\ = -9{k2MM)- Performing the Bo- 
goliubov transformation on H/^ is standard but lengthy. 
By manipulating indices under the sums we arrive at a 
relatively compact form for the part quartic in Bogoli- 
ubov operators: 



^ ^/ci+/c2+/c3+/c4,o|^o(^l, ^2, ^3, ^4) 



1,2,3,4 



4x42 43 44 + 



Vi{ki,k2, ks, /C4)4i<^-/C2«-/C3^-/C4 + h-c- + ^2(/^i, /^2, ks, /C4)4i 42^-^3^-/04}- (^^) 



The interaction vertices are given by 



Vo{ki,k2,ks,k4) = sgn(P) cos(A;p(i) - /cp(2) - ^fep(i) + ^/cp(2) + ^/cp(3) -^fep(4)), 

PeS4 



(22) 



with permutation P acting on the set {1, 2, 3,4}, 



Vi{ki, /C2, /C3, k^) =^y:^ Yl sgn(P)[sin(A:i - kp^2) - Ok, + Ok^^^^ + Ok^^^^ - Okp^^^) 
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■ sin(/cp(2) - kp^s) + ^/ci - ^/cp(2) + ^/cp(3) - ^/cp(4))] 



(23) 



with permutation P acting on the set {2,3,4} and finahy 

J . 



V2(ki,k2, fe, ^4) = - ( ^ sgn(P) cos(A;3 -k^- Ok^ + - Ok^^^^ + ^/cp(2)) 



PeSs 

+ ^ sgn(PO cos(/ci -k2- Ok, + Ok, - ^fc^,(3) + Ok^^^J 
P'eS3 

+ XI sgn(P)sgn(PO [cos(A;p(i) - kp^^s^ - Ok^^,^ + Ok^^,^ + 6>fe^,(3) - 6>fe^,(4)) 

PeS3 P'eS3 

where P and P' act on {1, 2} and {3,4} respectively. 



(24) 
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FIG. 1: (Color Online) The self-consistent Bogoliubov pa- 
rameter Ok, for J = 1. The parameter scales as A~^. 



New quadratic terms are generated by normal ordering 
the quartic piece. We must then include these with the 
original terms from H2 and solve for Ok self-consistently 
so that the off-diagonal terms are zero. This requirement 
may be recast as a self-consistency condition for every k: 



tan(2(9/e) 



2 sin(2fe) + f sin(fe) + ^ Q2{k, g) 
A + 2 cos(2fe) + f cos(^) + ^ Q) ' 



(25) 



where we have defined 



ei{k, q) = 2f^k,q)i-k-q) sin^(^g) + 9{k,q-q) sin(2(9g). 



(26) 



02(/^,g) = f{k-k){q-q) sin(26>g) - 4.g^k-k,q) sin^(^g). 



(27) 



Clearly the dispersion is also affected by the new 




FIG. 2: (Color Online) The single particle dispersion relation, 
efc (with J — 1) for A = 10. The dispersion calculated with 
and without self-consistency (solid and dashed curves respec- 
tively) and the exact spinon dispersion^^'^^ (dotted curves) 
are shown. The spinon result is not available in the case of 
a finite transverse field. For h = the dispersion is nearly 
sinusoidal. Note the functions are tt periodic for /i = and 
27V periodic otherwise. 



quadratic parts, becoming 

= j([^^ cos(2/c) + J cos(/c) ^^Yl ^) 

(28) 



h 



sin(2/c) + — sin( 



Evaluating the self-consistency and dispersion relations 
(25,28) numerically we can compare the gap (i.e. low- 
est excitation energy) to the mean field result found by 
Gomez-Santos.^^ Summing over 100 sites our results for 
the physical (two particle) gap are in excellent agreement. 
The self-consistent Bogoliubov parameter is plotted in 
Fig. 1 for a range of parameters. Figures 2 and 3 show 
that Ck is an excellent approximation to the spinon dis- 
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we write the commutator in the fermionic basis 



FIG. 3: (Color Online) The single particle dispersion relation, 
Ck (with J = 1) for A = 3. The dispersion calculated with 
and without self-consistency (solid and dashed curves respec- 
tively) and the exact spinon dispersion^^'^^ (dotted curves) 
are shown. 



persion. It is also apparent that use of a self-consistent 
Bogoliubov parameter is a small effect on the level of the 
dispersion, except in the presence of a transverse field. 

It is worth emphasising that the fermions that feature 
in the diagonalized quadratic Hamiltonian, Eq. (16), are 
not the same as the spinous of the exact treatment. ^^^^ 
Here fermion number is not conserved by the interac- 
tion vertices. In contrast spinon number is conserved by 
the exact solution. Instead the fermions described by 
should be viewed as propagating domain walls. 



3. Properties of the Eigenstates 

Previously it has been suggested that the fundamental 
excitations of Heisenberg-Ising chains are chiral.^^ We 
now make some remarks on this possibility, in light of 
our results. The relevant chiral operator, Cx^ is defined 
as 

= X • ^ 5„ X 5„+i = ^ SlSl^^ - S^^Sl^^. (29) 

n n 

For chirality to be a good quantum number for the XXZ 
chain, Cx must commute with the Hamiltonian. We find 

[H^ Cx] = ^ i{S^_i — iS^ + ^n+l ~ 



. ( qy qy 



(30) 



which is an 0(1) contribution. A priori this demonstrates 
that at least one eigenstate of the Hamiltonian is not an 
eigenstate of Cx- However it is possible to show that 
spinon states are generally not chirality eigenstates. First 



l,n+l 



l,n+l j 



4 ^ \-n-'^ 2 
n 

X (Mn-2,n-l - Mn-l,n + Mn,n+1 - Mn+l,n+2), (31) 

where we have defined 



^n,n+l = cl^li^l + CnCn+1 cJ,Cn+l CnC^+i- (32) 

As written, Eq. (31) contains only terms quartic and sex- 
tic in the creation and annihilation operators. However 
once Fourier transformed, rewritten in the Bogoliubov 
basis (a^) and normal ordered, quadratic terms will be 
generated. A true (one-spinon) excitation of the system, 
1^), involves a superposition of domain walls created by 
the al. operators. Schematically 

N 

l*) = E n h{h,--- ,h)A^-'al\0). (33) 

i—1 /ci ,ki 

where the /^(/ci,--- ^ki) are c-number functions of the 
momenta and we have written the small expansion pa- 
rameter A^~* explicitly. It has been shown that at low- 
est order the excitations are eigenstates of the chiral- 
ity operator. However from Eqs. (31) and (33) we see 
that the expectation (^| [H^ Cx] |^) will generally not be 
zero, instead having a finite contribution at lower orders 
in A~^. As A ^ oo these contributions vanish and at 
the Ising point the excitations can be chosen as chirality 
eigenstates. 



IV. DYNAMICAL RESPONSE 

The quantity of interest for inelastic neutron scattering 
is the dynamical structure factor S^^{uj, Q) given by 



(34) 



Here = \(7i is the a component of the spin operator 
at site /. The structure factor is related to the retarded 
dynamical susceptibility Xni^^Q) 



1 



1 



TT 1 - e-/5^ 



Im [x^^^,Q)] 



(35) 



r/5 



where 

JO 

x''\r,Q) = -1 ^e-^«('-'')(T.5r(r)5f,). (36) 



LI' 



Here we have introduced the Matsubara formalism 
and the expectation implies a thermal trace (• • • ) = 



7 



Following the discussion in Sec. II it is apparent that 
for = 0, x^^i^i Q) hence S^^{iJ, Q) will be diagonal 
in the indices a, b but that this is no longer the case for 

> 0, in agreement with experiment. 



A. Dynamical Structure Factor in the Fermionic 
Representation 



where 



^/3(^) =( §sin(7/e), cos(7/e), -§sin(7/e) )/3, (40) 
-fk=k^Q/2-0k-0k^Q, (41) 

and the 3x3 matrix n^^(r, Q\k^k') is given by 



As previously discussed, the Jordan-Wigner transfor- 
mation introduces non-local 'strings' which make calcu- 
lating complicated. We instead focus our attention 
on the transverse susceptibility, - 

First the required spin operator must be written in 
terms of the new fermionic operators: 

k 

-koik+q)]- 



- i sin(/c -Ok- Ok+Q + Q/2){ala 



-k-Q 



(37) 



We will use Eq. (37) to evaluate the time ordered dy- 
namical susceptibility. 



Up^{T,Q\k,q) = -{TrX^p{T,Q\k) Xt^(0,Q|g)[/(/3)), 

(42) 



^a-kak^Q 
Xf3j,{r,Q\k) = ( o^lak^Q 




aia"^ 



k^-k-Q/ 

(43) 

The imaginary time evolution operator in the interaction 
picture is 



U{r) 



--Trexp(^- £ dTiHi{n)y (44) 



X^^(T,g) = -^(T.cT5(r)a^Q). 



(38) The Fourier transform of the matrix 11 is given by 



As we aim to calculate (38) in perturbation theory in H4 
we now switch to the interaction picture. In order to 
simplify the perturbative calculation of x^^ it is useful 
to express (38) in terms of a 3 x 3 matrix 11/3^ (r, Q\k, k') 
(the matrix indices take values 7 = 1,2,3) as follows 



X^ 



(^'^) = ^E^/5(^)n/37(^.QI^.^0^:^(^0, (39) 



k,k' 



U{iUn,Q\k,q) = / dre'"^-^ Ii{T,Q\k,q). (45) 
B. Zeroth Order 



At zeroth order in perturbation theory we replace U {/3) 
in (42) by 1. All off-diagonal elements then vanish and 



we find {Unm = - ^m) 



Iil^{ii^Jn,Q\k,q) 
Iil2{^Un,Q\k,q) 
Iilr^{iUn,Q\k,q) 



\ X] Gidiii^nm, k + Q)GQ(iuJm, -k) [5k-q-Q - h,q] = ^ ^ ^fe+Q 1 [h-q-Q - , (46) 



lUn—ek—ek+Q 



\ Co(^^m, k + Q)Go{-iLUnm, k)Sk,q = .^^^ ^^^^^ „ ^k,q, 



ILdn^^k — ^k+Q 



(47) 



- ^ Go(-itJnm, -k - Q)Go{-iu;m, k) [5k-q-Q - 5k,q] 



rik + ^fe+Q - 1 



[^k,q — ^k,-q-Q] • 



(48) 



Here Uk = l/{e^^'' + 1) and the bare Green's function is 
given by 

Go{ikn,k) = — ^ . (49) 

IKfi Ck 

The dynamical susceptibility at zeroth order in perturba- 
tion theory is then obtained by substituting the matrix 



into (39) and carrying out the momentum sums. Tak- 
ing the thermodynamic limit and analytically continuing 
to real frequencies, icOn uo^ir]^ we arrive at the follow- 
ing expression for the zeroth order retarded susceptibility 
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^ dk 

8^ 



{1 - cos{2k ^ Q - 2[0k ^ Ok+Q])) 



2(1 + cos(2/c + Q - 2(^fc + ^fc+g))) 



(50) 



The remaining /c-integral cannot be carried out analyt- 
ically in general as the Bogoliubov parameters Ok need 
to be determined self-consistently and are therefore only 
known implicitly. However, in the limit A ^ oo the inte- 
gral can be taken and simple expressions for Xrq){^->Q) 
may be obtained. 



1. The A ^ oo,T = Limit 

In order to evaluate the susceptibility further we take 
h = and expand Ck as a series in 1/A. This gives 



Ck = ^ — h Jcos(2/c) + sin(2/c) 



J(A + 2 cos(2fc + Q) cos(Q)) 
efc+Q - Efe = -2 J sin(2A; + Q) sin(Q) + . . . 



We see that poles occur at 



LO = 2Jsin(2A;o + Q) sin(Q), 

uj = J(A + 2 cos(2A;_ + Q) cos(Q)), 

w = - J(A + 2 cos(2fc+ + Q) cos((5)). 



(51) 

...,(52) 
(53) 



(54) 
(55) 
(56) 



The contribution from Ok only enters at O(^) so we ne- 
glect it. We wish to take the imaginary part of the re- 
tarded susceptibility as this is proportional to the dy- 
namical structure factor. Defining 



P{u,E) = Q{\E\^io)e{\E\-io) 



if |a;|>|^| 
if Iwl < \E\ 



(57) 



where the 6's are Heaviside step functions, we find 



(l + cos(2/co + g)) 
1 



(^feo+g - nko)P{uj, 2J sin(Q)) 
16J|cos(2A;o + Q)sin(Q)| 

{rik+^Q + rik+ - 1) 



^(l-cos(2/c+ + Q))_^, . 
2^ ^ ^ ^^^16J|sin(2A:+ 

X P(ct; + AJ,2Jcos(Q)) 



-g)cos(Q)| 



-(1 -cos(2A:. 



■Q)h 



{nk_+Q +^/e_ - 1) 



16 J| sin(2fe_ 
X P(cj- AJ,2Jcos(Q)), 



g)cos(Q)| 

(58) 



with 



2/co + Q = arcsin 
2/c+ + Q = arccos 
2k- + Q = arccos 



2Jsin(Q) 
1 



2cos(Q) \J 
1 



_2cos(Q) Vj 



(59) 
(60) 
(61) 



Using the expansion of e/e one can also show that at next 
to leading order 



sinh(^/3cj)P(cj,2Jsin(Q)) 



cosh(^^cc;) + cosh(f JA - f cot(Q)((2Jsin(Q))2 - cj^y/^) ' 



cosh{^ poo) P{u; ± JA,2Jcos(Q)) 
cosh(^/3cj) +cosh(f tan(Q)((2Jcos(Q))2 - {cj ± JA)2)i/2) * 



nk+-\-o + ^/e+ — 1 = ^ ^ ; — • (62) 
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Inserting these relations into the fuh expression yields 



1 



1 



((2Jsin(Q))2 -cj2)i/2 2Jsin(Q) 
1 /2Jcos(Q) + (ct; + JA) 



sinh(^/?cj)P(cj,2Jsin(Q)) 



cosh(^/3cj) +cosh(fJA 



f cot(Q)((2Jsin(Q))2 



.)l/2) 



cosh{^ poo) P{uj + JA, Jcos(Q)) 



32Jcos(g) V2Jcos(Q) - + JA)y cosh(^/3cj) + cosh(f tan(Q)((2J cos(Q))2 - (tj + JA)2)i/2) 

1 /2Jcos(Q) - (cj - JA)y cosh(^^cj)P(tj - JA,2Jcos(g)) 

32Jcos(Q) V2Jcos(g) + (cj- JA); cosh(^/3c^) + cosh(f tan(g)((2J cos(g))2 - {u - JA)2)i/2)' 



(63) 



This anaysis shows that, at zeroth order, the response 
diverges as an inverse square root both for uo = ±(AJ — 
2Jcos(g)) and uo = ±2J sm{Q). These divergences are 
associated with the gapped two-spinon response and the 
thermally activated Villain mode respectively. 



C. First Order Perturbation Theory 

At first order two classes of contributions appear, 
which may be distinguished by considering their diagram- 
matic representation. As we have already seen, the three 
zeroth order diagrams take the form of bubbles. The 
first class of diagrams in first order consists of two bub- 
bles connected by an interaction vertex. The second class 
of diagrams consists of a single bubble with a self-energy 
insertion. We detail these contributions in Appendix B. 
An important feature is that several of the first order 
contributions have stronger singularities as functions of 
the external frequency and momentum than the zeroth 
order results. This indicates that it is necessary for the 
expansion to be resummed before useful physical results 
can be extracted. 



V. BUBBLE SUMMATION 

We have shown that at zeroth order the susceptibility 
diverges for certain co and Q and that this is matched 
by stronger divergences in some of the first order terms 
(see Appendix B). Moreover, it is clear from the first 
order calculation that higher orders in perturbation the- 
ory will exhibit stronger and stronger divergences. In 
order to get physically meaningful results we therefore 
should sum the most divergent classes of diagrams. In 
the case at hand, the complicated momenta dependence 
of the vertices and the self-consistent Bogoliubov trans- 
formation makes the determination of the most divergent 
contributions an impossible task. Instead we resum just 
the connected bubble diagrams and justify our choice by 
comparing results with the exact T = calculation. We 
explain how to incorporate self-energy corrections in Sec. 
VD. 



A. RPA-Like Scheme 

The RPA scheme consists of carrying out bubble sums 
of the type shown in Fig. 4 





+ ••• 



FIG. 4: Bubble summation for one contribution to the dy- 
namical susceptibility matrix x^'^(*<^n, Q)- The thick lines in- 
dicate that the single particle propagators may be resummed 
to include self-energy corrections. 



To carry out these summations is non-trivial as the lines 
in the internal bubbles can take any orientation and the 
momentum dependence of the vertices is very compli- 
cated. In order to proceed we organize the interaction 
vertices into a 3 x 3 matrix 

Vn{Q\k,q)=V2ik + Q,-k,q,-q-Q), 
Vi2{Q\k, q) = -Wi{q + Q, -q, k, -k - Q), 
Vi3{Q\k, q) = 6Vo{k + Q, -k, q, -q - Q), 
V2i{Q\k,q) = 3Vi{k + Q,-k,q,-q - Q), 
V22{Q\k, q) = -AV2{k + Q, q, -k, -q - Q), 
V23iQ\k, q) = -Wiik, q + Q, -q, -k - Q), 
Vzi{Q\k,q) = Wo{q + Q,-q,k,-k - Q), 
V32{Q\k,q) = 3Vi{q,k + Q,-k,-q-Q), 
V33{Q\k, q) = V2{q, -q-Q,k + Q, -k). (64) 
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As is shown in Appendix C, in the thermodynamic hmit 
the RPA-hke bubble summation without taking into ac- 



J 



count self-energy corrections results in an integral equa- 
tion of the form 



n^^A(icu„, Q\k, k') = n%{iun, Q\k, k') 



J Jx„^(iw„, Q\k, q)n^^^{iu;n, Q\q, k') . (65) 



where the kernel K is defined as the infinite volume limit Defining a convolution * by 
of 

KapiiiOn, Q\k, q) = ^n°^(iw„, Q\k, k')V-,p{Q\k' , q). 

k' 

(66) 



{X*Y) 



(67) 



r 



this can be rewritten as 

(/ - K) * n 

The integral equation (68) is then readily solved 

nf^PA = (j_K)"'*no. 

After analytic continuation to real frequencies 
irj the quantity of interest is calculated as 

dk dq . 



RPA 



(68) 



(69) 



. (2^)^ 



-L^{k)nll^{LO + iv, Q\k, q)L\{q). (70) 



In practice the solution of the integral equation discussed 
above is reduced to a simple matrix inversion problem. 
We discretize all momentum integrals in terms of sums 
over N = 400 points. We find that this value is large 
enough to make discretization effects negligible. We set 
J = 1 and the regulator 77 = 10~^. The discretized rep- 
resentation of J — * V (which is a matrix both in 
momentum space as well as in the 3x3 matrix space 
labelled by greek indices) may be found using standard 
linear algebra routines. The dynamical structure factor 
is then evaluated at 24000 points in frequency space. Us- 
ing a finite system size results in a finite number of poles 
in the susceptibility (36). In turn, because r] is finite, the 
calculated structure factor will be composed of a number 
~ of Lorentzian peaks of width 77. Finally we convolve 
this result with a suitable Gaussian. 



B. Comparison with Exact Results for T — 

Given the uncontrolled nature of our bubble summa- 
tion it is essential to compare it to exact results at zero 



: 0.04 



1 1 1 1 
Q=0, A=10, T=0 /v '' 


1 ' 1 ' 1 ' 

— Resummed 
---- Exact 
%N IS 


■iT 1 1 1 


III, r-\ , 



FIG. 5: (Color Online) The dynamical structure factor as 
found by resummation, the exact result^^'^^ and the calcu- 
lation of Ishimura and Shiba^° (IS) at T = 0, Q = and 
A = 10. In all cases the curves are convolved with a Gaus- 
sian in frequency space of full width half maximum 0.12. 



temperature^^ '^^ in order to assess its quality. In Figs. 
5, 6, 7 and 8 we plot our results against the exact results 
for the dynamical structure factor for A = 10, T = 
and several momenta. We also include for comparison an 
earlier result for the DSF due to Ishimura and Shiba.^^ 
We see that for A = 10, T = and at most wave 
vectors the resummation is a highly accurate approxima- 
tion. This suggests that the diagrams we have selected 
account for most of the spectral weight at low temper- 



11 





1 ' 1 


1 , 1 , 
— Resummed 


Q=7c/4, A=10, T=0 // ^ 


<^ ---- Exact 






\ IS 
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i 1 , 1 


1 \ \ 1 
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CO 

FIG. 6: (Color Online) The dynamical structure factor as 
found by resummation, the exact result ^^'^^ and the calcu- 
lation of Ishimura and Shiba^*^ (IS) atT = 0, (5 = 7r/4 and 
A = 10. In all cases the curves are convolved with a Gaussian 
in frequency space of full width half maximum 0.12. 




FIG. 8: (Color Online) A comparison of our calculation and 
the exact result^^'^^ for Q = 7r/2 at T = and A = 10. Left 
panel: 7r/2 curves convolved with a Gaussian (width 0.08) 
plotted with the exact result at Q = 0, in order to demon- 
strate the scale. Right panel: The same tt/2 curves but plot- 
ted at a different scale. At Q — 7t/2 the result of Ishimura 
and Shiba is undefined as the denominator goes to zero.^° 



0.12 



0.1 



0.04 



0.02 
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— Resummed 


Q=7C, A=10, T=0 


\ ---- Exact 
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FIG. 7: (Color Online) The dynamical structure factor as 
found by resummation, the exact result ^^'^^ and the calcu- 
lation of Ishimura and Shiba^° (IS) at T = 0, Q = tt and 
A = 10. In all cases the curves are convolved with a Gaus- 
sian in frequency space of full width half maximum 0.12. 



ature. It is known from the exact result^^'^^ that the 
gapped response diverges along its lower energy thresh- 
old, in a region of momentum centred about 7r/2 (the 
size of which increases as A 1). Correspondingly the 
resummation is less accurate near 7r/2, although there is 
still qualitative agreement (Fig. 8). For larger values of 
A our approximation becomes even better. 



C. Analytical Resummation for T = 

It is possible and instructive to evaluate the resumma- 
tion exactly in the limit T = 0, A — > oo. At zeroth order, 
positive frequencies and T = 0, the only non-zero dia- 
gram is the particle-particle propagation bubble. This 
diagram leads to a response in the region co ^ AJ. As 
T ^ the thermal occupation factors vanish and by ex- 
panding in 1/A we find 



(71) 



so that we can neglect the Bogoliubov phases. The dia- 
gram then reduces to 



dk sin^(A: + Q/2) 



A{iu;ri,Q)' 



Att iun - e/c - e/e+Q 
Using the approximation 

e/e + e/c+Q ^ A J + 2Jcos(Q) cos(2/c + Q) 



(72) 



we have 



^ dk sin2(/c + g/2) 
1 



/' 

Jo 



dk 1 — cos(2/c - 



Q) 



-AJ 



4Jcos((5) Jo 27r a; - cos(2/c + Q) 
The remaining integral can be taken 



and u = 75-7 — 77^. 

2 J cos(Q) 

by standard contour integration methods and analytic 
continuation to real frequencies is then straightforward. 
With the approximations made here the vertex V2 takes 
a particularly simple form and resummation amounts to 
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the geometric series 

-Im[x""(w, Q)] = - lm{A{w, Q) + 4Jcos(Q)^(a;, Qf 
+ (4J)2cos2(Q)A(w,Q)3 + ...) 
A{lo,Q) 



Finally the T = result, to lowest order in A is 



Im( 



1 -4Jcos(Q)A(cj,Q) 



). (73) 



r V^(2Jcos(g))2-(a;-AJ)2 





AJ| < |2Jcos(Q)| 
otherwise. 



(74) 



r 



Ishimura and Shiba calculated the DSF at T = using 
a method based on perturbation theory combined with a 
Lehmann representation.^^ The expression obtained here 
(74) agrees to lowest order with their result. 



D. One— Loop Self— Energy Corrections 

As they stand, the first order tadpole contributions, 
Eqs. (B10-B16), cannot be incorporated into the RPA. 
This is because they contain divergences which are not 
resummed according to the scheme above. Instead we 
must first calculate a new single particle Green's func- 
tion, using a Dyson equation to resum the tadpole self- 
energy corrections. This Green's function is then used to 
calculate a new bubble diagram. 

Including the tadpole corrections generates anomalous 
propagators at higher orders. These are most efficiently 
taken into account by formulating the propagator as a 
matrix. 



g{iujn,k) 



-f 

Jo 



dTe"^"^giT, k) , 



At zeroth order the result is 







(75) 
(76) 



The single particle self-energy is defined by the Dyson 
equation 

g~^{iun,k) = gQ^{iUn,k) - T>{iUn,k). (77) 

If only 'tadpole' type diagrams are included the self- 
energy is, to first order in perturbation theory, frequency 
independent 

^ ^ 77^ / W2{k,p, -p, -k) -3Vi(p, -p, k, -k) 
2^ N\ 3Fi(p, -p, k, -k) -W2{k,p, -p, -k) 

(78) 



The elements of the matrix are (9(n/cA^) and hence their 
effects will be most pronounced when the gap is small or 
the temperature is large. 
Using the relations 

gl^{iujn^k) 
gl^{iun,k) 

si 

one finds that 
giiujn.k) 



1 




1 

- e/c 


1 






1 


(E-) 









(79) 

(80) 

(81) 
(82) 



1 



(icu„)2-(efc+Sii)2-|Si2|2 

11 _5^12 

iuj„ - ek- 



X 



lOJr, 



h efe + Si 
Si2 



^11 



We rewrite this as 



^^(iujn.k) 



ILOn — Ek 

Z7 



^^{iun.k) 



• Ek 
1 



lUJri 



■ Ek- 
1 



ILdn + Ek ILdn 
12/ 



Ek 



g (iujn,k) = -g {iun.k), 
with the definitions 

i?fe = V(efe + S")2 + |Si2|2, 



It 



efc 



Ek 



Si2 
2Ek 



(83) 
(84) 

(85) 
(86) 

(87) 

(88) 



E. 



Bubble Summation with Self— Energy 
Corrections 



The one-loop self-energy corrections to the propaga- 
tors can be taken into acount in the bubble summation 
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for the dynamical susceptibility as follows. We define a 
3x3 matrix n*^ by 



nl^{T,Q\k,q) = 

-{T,Xf,p{T,Q\k) X\^{0,Q\q)U{/3)) 



1-loop S 



,(89) 



where only one-loop self-energy corrections are taken 
into account. This amounts to calculating the two-point 



function of X and using the (anomalous) propagators 
(82). The elements of n*^ are listed in Appendix D. We 
now follow the same steps as in Appendix C and show 
that summing all bubble diagrams of the form shown in 
Fig. 4 gives rise to an integral equation obtained from (65) 
by the replacement E^^ 



I.e. 



J 



U^$^{^uOr^. Q\k, k') = nf^(ic^,, Q\k, k') 



I 



dq 



2^-„^(ia;„, Q\k, q)n^J^(iw„, Q\q, k') . 



(90) 
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FIG. 9: (Color Online) The dispersion Ek^ including tadpole 
self-energy corrections for A = 10. The left panel shows 
the gradual narrowing of the bandwidth with temperature. 
At this scale the dispersion for T = J is indistinguishable 
from that for T = 0. The right panel shows the two distinct 
elements of the self energy matrix, Xlgp at T = 2 J. 



where the kernel is defined as the infinite volume 
limit of 

K^piiLOn. Q\k, q) = J2 nf^(^^"' Ql^' k')V^0{Q\k', q). 

k' 

(91) 

At T = the one-loop corrections to S) vanish, so that 
we recover our previous result. On the other hand, for 
T > 0, X) plays an important role, altering the dispersion 
and shifting the thresholds of the dynamical response (see 
Fig. 9). 

One can also calculate the two-loop corrections to the 
single particle propagator but it is not as simple to in- 
corporate them into the RPA. For reference they are in- 
cluded in Appendix E. 



VI. RESULTS AND DISCUSSION 

We now turn to a discussion of our results at finite tem- 
peratures and applied fields. For the perturbative expan- 
sion to be valid we require A 1. In selecting diagrams 
to resum we have favoured those that are most relevant 
at low temperatures, i.e. those that feature few thermal 
occupation factors. The relevant energy scale is the sin- 
gle particle gap ~ A J/2, so we restrict our discussion to 
temperatures T < A J/2. We calculate the transverse dy- 
namical structure factor as described in Sec. VA, using 
the matrix XI^ found in Sec. VD and setting A = 10, 
J = 1. For A 1 the dynamical structure factor (at 
positive frequency) consists of two pieces: a gapped con- 
tinuum response occuring at frequencies uo ^ A and a 
response for a; ~ that is only seen at finite tempera- 
ture. 

On general grounds one expects that at finite temper- 
ature the very sharp thresholds seen at T = should dis- 
appear. In our approach the thresholds are still present 
although they are obfuscated by the necessity of convolv- 
ing the response with a Gaussian. This is a consequence 
of the diagrams we have taken into account. We expect 
this 'thermal broadening' to be a small effect that could 
be taken into account by including certain two-loop di- 
agrams. These two-loop diagrams connect the response 
to decay channels of higher particle number. We discuss 
the issue further in Appendix E. 



A. Gapped Response for /i = 

We first consider the gapped (interband) response at fi- 
nite temperatures and for = 0. At T = the transverse 
DSF is dominated by a two-spinon scattering continuum 
that occurs at energies around twice the single-particle 
gap, i.e. uj ~ A J. In Figs. 10, 11 and 12 we show how 
the DSF in this regime of energies changes at finite tem- 
perature. The most striking feature is a narrowing of the 
response with increasing temperature. This is in agree- 
ment with inelastic neutron scattering experiments on 
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FIG. 10: (Color Online) The dynamical structure factor at 
finite temperature for a; A and A = 10 at wave vector 
Q = 0. 



FIG. 12: (Color Online) The dynamical structure factor at 
finite temperature for cj ~ A and A = 10 at wave vectors 
Q = 7r/4 and 7r/2. Note that the scale of the intensity axes 
differs dramatically between the two plots. 




CO 



FIG. 11: (Color Online) The dynamical structure factor at 
finite temperature for a; A and A = 10 at wave vector 

Q = 7T. 



TlCoCls^^. Another notable feature at T = is that the 
response is not symmetric about Q = 7r/2.^^ As can be 
seen in Figs. 5-7 our calculation captures this behaviour. 
Accordingly the response develops asymmetrically with 
temperature, and in a non-trivial way. In particular the 
maximum of the response for < Q < 7t/2 moves to 
lower frequencies as temperature increases, but is shifted 
to higher frequencies for 7r/2 < Q < tt. 

For all wave vectors the total spectral weight in the 
gapped region decreases as temperature increases. The 
thresholds of the response shift as temperature increases 
due to the thermal dependence of the single particle dis- 
persion Ek, Eq. (87). This depletion of spectral weight is 
physically sensible because the excitations are fermionic; 
as temperature increases more states are thermally oc- 
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FIG. 13: (Color Online) The dynamical structure factor for 
a; and A = 10. Main plot: the Villain mode for a range 
of wave vectors. Inset: the development of the Villain mode 
at Q = TT with temperature. 



cupied and concomitantly there are fewer states for the 
new pair of fermions to fill. 



B. Villain Mode forh = 

At temperatures greater than zero there is a thermal 
population of spin excitations. Neutron scattering can 
then lead to processes that do not change the spinon 
number of a given microstate, thus giving rise to an in- 
traband respose at low energies, a; ~ 0. Dynamics of 
this kind were first described by Villain for the case of 
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FIG. 14: (Color Online) Comparison of the T = dynami- 
cal structure factor for A = 10 with and without an appUed 
transverse field h = J/2. Main panel: the divergent response 
at 7v/2 is split into a much broader two peaked structure by 
the application of the field. Inset: the effect of the field at 
Q = is only apparent at the upper threshold. 



an finite length chain with an odd number of sites. As 
such processes rely on states being thermally occupied, 
their contribution to the DSF grows with temperature. 
The principle feature of the response is a well defined 
resonance or 'mode' at a; = 2Jsin(Q).^'^'^^ At lowest or- 
der in our calculation the intraband response exhibits a 
square-root divergence (IV B 1). Taking interactions into 
account by our bubble summation leads to a smoothing 
of the divergence, which however occurs in a very small 
region in energy, close to the threshold. This means that 
the zeroth order (especially when convolved with an ex- 
perimental resolution) is an excellent approximation to 
the resummed result. 

The intraband response is also asymmetric about Q = 
7r/2 (see Fig. 13). For < Q < 7r/2 the response between 
the peaks is suppressed relative to that at Q = 7r/2. For 
7r/2<(5<7ritis enhanced. 



FIG. 15: (Color Online) Wave vector dependence of the dy- 
namical structure factor for A = 10 at /i = J/2, T = 0. The 
main panel shows wave vectors Q = 0, 7r/4 and 7r/2. The inset 
shows plots for Q = 37r/4 and tt. 
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FIG. 16: (Color Online) Wave vector dependence of the Vil- 
lain mode for A = 10 at /i = J/2,T = J. The main panel 
shows wave vectors Q = 7r/4,7r/2 and 37r/4. The inset shows 
the response at Q = 7r/2 for /i = and /i = J/2. 



C. Response in a Transverse Field 

The transverse field h only enters the quadratic part of 
the Hamiltonian (16), hence its influence on the scatter- 
ing response is through the single particle dispersion and 
not the interaction vertices. The first property to take 
into consideration is that h will have an effect on the ex- 
citation gap. As the magnitude of the field h approaches 
AJ the gap collapses and the perturbative expansion is 
inapplicable. Hence we consider a field small compared 
to A J. For fields A J the gap opens again with the 
chain ferromagnetically ordered. A second important fea- 
ture to note is that a non-zero h causes the period of 
to double (see Fig. 1): the maximum at /c = (2n + l)7r is 



reduced relative to those at /c = 2n7r. This leads to the 
double peaked structure seen at Q = 7r/2 in the gapped 
response (Figs. 14, 15) and throughout the low energy 
scattering (Fig. 16). This splitting of the Villain mode 
peaks was observed by Braun et al}^ 

The asymmetry of the response about Q = 7r/2 is 
further increased by the transverse field. In particular, 
though the exact result at T = Q^^'^^ and A = 10 shows 
that the energy thresholds of the gapped response are 
very nearly symmetric about 7r/2, this is no longer the 
case in a finite field. Instead the upper threshold near 
Q = is pushed to higher energies but is relatively un- 
changed near Q = tt. We also note that the narrowing in 
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energy of the response near Q = 7r/2 is suppressed rela- 
tive to h = 0. Asymmetry is also seen in the thresholds 
of the low energy response (Fig. 16). 



VII. COMPARISON TO DIAGONALIZATION 
OF SHORT CHAINS AT T > 

As we have seen above, at zero temperature our ap- 
proach gives good agreement with the exact DSF. In or- 
der to assess the quality of our approximate DSF at finite 
temperatures, we have computed the DSF by means of 
numerical diagonaliztion of the Hamiltonian on finite pe- 
riodic chains of up to 16 sites. The dynamical suscepti- 
bility (36) is expressed by means of a Lehmann expansion 
in terms of Hamiltonian eigenstates |n) of energy En as 

n,m 

X (n| 1"^) ("^1 ^0 In) , (92) 

where the sums are taken over the eigenbasis of H. For 
a sufficiently small system the eigenstates can be calcu- 
lated numerically. Due to the exponential increase in 
the dimension of the Hilbert space with the number of 
sites the method is restricted to short chains. We con- 
sider systems with N < 16. In order to approximate 
the thermodynamic spectrum from a finite number of al- 
lowed transitions, we take the parameter r] in (92) to be 
sufficiently large, so that the finite sum of Lorentzians 
in (92) becomes a smooth function of (the real part of 
the) frequency. Clearly this procedure can give a mean- 
ingful approximation of the susceptibility in the thermo- 
dynamic limit only if r] is very small in comparison to 
the thermal broadening. Hence the method is restricted 
to sufficiently large temperatures. In order to obtain a 
measure of the importance of finite-size effects we calcu- 
late the DSF for system sizes = 8, 12, 16. We find that 
the finite-size effects depend strongly on which region of 
energy and the momentum we consider. 

In Figs. 17-20 we compare the dynamical structure 
factor calculated for a 16-site chain to the results ob- 
tained by our perturbative approach for several values of 
momentum. For Q = 0, tt and T = 2J (Fig. 17) the 
agreement of the two methods is good. The difference 
at very small freqencies is probably due to finite-size ef- 
fects in the exact diagonalization results. The finite-size 
effects for the main peak are found to be quite small. 

At Q = 7r/2 and T = 2J (Fig. 18) the agreement of 
the two methods is less impressive. The disagreement at 
high frequencies is likely due to inaccuracy of the bubble 
summation in our perturbative method (we recall that 
the agreement of our resummation with the exact result 
at T = was worst for Q = 7r/2). On the other hand, 
the exact diagonalization results are found to still suffer 
from finite-size effects at small frequencies. 
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FIG. 17: (Color Online) Comparison of the dynamical struc- 
ture factor obtained from the perturbative approach to exact 
diagonalization on a 16-site chain with r/ — 0.5 J for T — 2J 
at momentum Q = 0. In order to facilitate a comparison the 
perturbative result has been convolved with a Lorentzian of 
width 77. 
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FIG. 18: (Color Online) Comparison of the dynamical struc- 
ture factor obtained from the perturbative approach to exact 
diagonalization on a 16-site chain with 77 = 0.5 J for T = 2J 
at momentum Q = 71. In order to facilitate a comparison the 
perturbative result has been convolved with a Lorentzian of 
width 77. 



VIII. CONCLUSIONS 

We have calculated the transverse dynamical structure 
factor of the XXZ spin chain at finite temperature and 
applied field. The perturbative method we use is ac- 
curate at low temperatures and for large A > 10. In 
this case the chain is in the Ising phase and the excita- 
tions are descended from propagating domain walls. The 
scattering response is composed of two distinct parts for 
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FIG. 19: (Color Online) Comparison of the dynamical struc- 
ture factor obtained from the perturbative approach to exact 
diagonalization on a 16-site chain for T = 2 J at momentum 
Q = 7r/2. Left panel: low frequency response with 77 = 0.2 J. 
Right panel: high frequency response with 77 = 0.05 J. In or- 
der to facilitate a comparison the perturbative results have 
been convolved with a Lorentzian of width 77. 




2. The position of the peak (as a function of fre- 
quency) of the high frequency (gapped) continuum 
at T > becomes asymmetric about Q = 7r/2 as 
temperature increases. 

3. The thermally activated response at low frequencies 

~ is asymmetric about Q = 7r/2. 

4. The Villain mode splits into two peaks in a trans- 
verse field. 

The main advantage of our method compared to previ- 
ous theoretical approaches lies in the fact that it is not 
restricted to asymptotically large values of A and treats 
the Villain mode and the gapped response in a unified 
way, which allows us to determine the ratio of spectral 
weights between these two features. Our results are in 
qualitative agreement with inelastic neutron scattering 



experiments. 



8,9,11,13 



It would be interesting to perform 



quantitative comparisons with polarized neutron scatter- 
ing data. 
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APPENDIX A: DIRECT JORDAN-WIGNER 
TRANSFORMATION OF THE HAMILTONIAN 

We may write the Hamiltonian (1) as H = Hq + i^i, 
where 



FIG. 20: (Color Online) Comparison of the dynamical struc- 
ture factor obtained from the perturbative approach to exact 
diagonalization on a 16-site chain with 77 = 0.05 J for T — 10 J 
at momentum Q = 0. In order to facilitate a comparison the 
perturbative result has been convolved with a Lorentzian of 
width 77. 



A > 10, namely a gapped continuum at ~ A J and 
the 'Villain mode' at cc; ~ 0. Our results pertain to the 
low temperature and field dependence of both. Our main 
observations are 

1. The response associated with the gapped two- 
spinon continuum at T = narrows in energy as 
temperature is increased and loses spectral weight 
to the emerging low frequency Villain mode. 



(Al) 



Hq can then be expressed as a quadratic form in spinless 
fermions by means of the Jordan- Wigner transformation 



1 , 



y 



(A2) 



The full Hamiltonian takes the form 
J 

4 



lE(A-i) 

n 



cUl+i + h.c. 



(A3) 
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The terms quadratic in fermions can be diagonalized by 
a Bogoliubov transformation (15). The resulting free- 
fermion dispersion is different from (16). Expressing the 
interaction part in terms of Bogohubov fermions and nor- 
mal ordering self-consistently leads to a theory of the 
same structure as the one derived in section III. In fact 
we expect it to be identical, but we have not verified this. 



APPENDIX B: ONE VERTEX DIAGRAMS 

In this section we list the first order contributions to 
the transverse response and discuss their divergent be- 
haviour. 



1. 'Connected Bubble' diagrams 



The contributions to L^U^^^cj, Q)R'^. given by connect- 
ing two bubbles are (the box vertices in the diagrams in- 
dicate the inclusion of the external factors, and R^): 




^2 ^^^(^fe) sm(7g)^2(-^, k^Q,q,-q- Q) 



7V2 



k,q 



Uk + n/c+Q - 1 ng + Uqj^Q - 1 

iuJn — ^k — ^k-\-Q 'l^n — ^q — ^q+Q 



(Bl) 




^ sin(7fe) sin(7g)l/o(-^, k^Q,q,-q- Q) 



k,q 



Uk + n/c+Q - 1 ng + riq^Q - 1 

icJn+e/e+efc+Q i^n — ^q — ^q+Q 



(B4) 




= jV2 S ^^'^i'yk) sin(7,)Vb(-fc, k + Q,q,-q- Q) 



k,q 



Uk + nfc+Q - 1 ng + Uq^Q - 1 



(B5) 




= ^ ^^^(^^) sin(7g)iFi(A;, g + Q, -q, -k - Q) 



k,q 



Uk - UkJrQ ^q + ^g+Q - 1 



(B6) 





51 sin(7fe) sm{-fq)V-2{q, -q - Q, -k, k + Q) 



k,q 



Uk + nfc+Q - 1 ng + Uq^Q - 1 

ZCJn + e/e + e/e+Q ZC^n + Cg + eg+Q 



(B2) 




^ sin(7/e) cos{-fq)iVi (-g, q^Q,-k-Q, k) 



k,q 

Uk + n/c+Q - 1 Uq- Uq^Q 
iuJn — ^k-^k+Q i^n — ^q-^^q+Q 

-q 




(B7) 



cos(7g)l/2(/^ + Q, q, -k, -q - Q) 



k,q 



rik — rik+Q 



^q+Q 



iuJn + ^k-^k+Q iuJn^eq-Cq^Q 



(B3) 



^ cos(7fc) sin(7,)iVi*(fc + Q, -k, q, -q - Q) 



k,q 



rik - rik+Q Uq + ng+Q - 1 

iuJn + ^k-^k+Q i^n-^q-^q+Q 



(B8) 
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k,q 

rik + Uk+Q - 1 Uq- Uq^Q 



(B9) 



2. 'Tadpole' type diagrams 

These contributions consist of bubbles in which one 
of the propagators features a tadpole type self energy 
interaction: 





= iV2E(^°«(27fe)-l) 



k,q 

Uk + n/e+Q - 1 

JuJn — e/c — e/c+Q 



nqV2{-k,q,-q,k) 
iiUn — e/c — e/c+Q 



(BIO) 





2 



2^(cos(27fc) - 1)^-— 



k,q 

Uk + n/e+Q - 1 



ILJn + e/e + e/c+Q 



(Bll) 





= —2 2^ COS 27fc + 1)^ — 



+ ^n/e(l - Uk) 



ILOn + e/e — e/e+Q 

^ Y^/ /o \ , .\'^q^'2.{k^Q,q,-q,-k-Q) 
— 77 > (cos(27^) + 1)— ^ 



k,q 



rik — ^/c+Q 



zcjn + e/c — e/c+Q 



(B12) 




7V2 



^ sin(27/e)zl/i (g, -g, -/c, k) 



k,q 



ILdn^^k — ^k+Q 



rik + n/e+Q - 1 2n(efe) - 1 



_ lUJn — ^k—^k+Q 



2e. 



(B13) 




3 

— ^ sin(27fc)iFi (g, -g, 

rik - rik+Q 2n(efc) - 1' 



Jp X! sm(27fc)^ Vi (g, -q, -/c, k) 



(B14) 



k,q 



lUJn — ^k—^k+Q 



rik - rik+Q 2n(e/e) - 1 



iuJn^^k — ^k+Q 



2e. 



(B15) 




^ ~ 7V2 ^ sin(27/e)i Vi (g, -q, -k, k) 



k,q 



iuJn — ^k^^k+Q 



rik + ^fe+Q - 1 2n(efe) - 1 



2efe 



(B16) 



One expects that, just as at zeroth order, these diagrams 
will have divergences at certain energies. Inspection of 
the pole structure of the connected bubble diagrams (Bl- 
B9) suggests that they will have stronger divergences 
than the zeroth order because they feature products of 
poles. Numerical results support this assumption for Eqs. 
(B1,B2). On the other hand, it is clear that Eqs. (B4- 
B9) will be very small (for large A) because they contain 
products of terms that, individually, are only significant 
for different discrete regions in uo. The behaviour of Eq. 
(B3) is subtler. This diagram gives the most significant 
first order contribution to the response at ~ 0. How- 
ever it does not contain a stronger divergence than its 
zeroth order equivalent. The reason is as follows: con- 
sider the sum 



E: 



I{k,Q) 



■ efc - Cfe+Q + iri 



(B17) 



with /(/c, Q) an analytic function of /c, Q. As shown at 
zeroth order, because the dispersion e/e is bounded, the 
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imaginary part of the sum (as 77 ^ 0) will in turn be 
bounded. In general there is a divergence at the thresh- 
old uo = max(e/e — e/e+g) of the corresponding response. 
Naively one then expects that in the double sum 

I{k,Q)I{q,Q)I\k,q,Q) 
^ (cj + e/e - e/e+Q + ir]){{uj + e/e - e/e+Q + ir])) 

with r{k^ Q) analytic, the maximum contribution will 
occur for A: = g and uj = max(e/e — e/e+g). For Eq. (B3) 
the function V2{k + Q^k^—k^—k — Q) vanishes at the 
threshold uo = max(e/e — e/c+g). This means that the di- 
vergence is in fact substantially weaker than at zeroth 
order. This behaviour does not persist when self-energy 
corrections to the propagator are included. This is be- 
cause the self-energy corrections shift the thresholds of 
the response. Stronger than leading order divergences 
are also found in Eqs. (B10-B16). To take them into ac- 
count the single particle propagator must be resummed 
using a Dyson equation, as in Sec. VD. 

APPENDIX C: MATRIX STRUCTURE OF THE 
BUBBLE SUMMATION 

In this appendix we show that summing all bubble di- 
agrams results in the integral equation (65) for the ma- 



trix n^^^. The proof follows by induction. Our starting 
point is expressions (39) and (42) for the dynamical sus- 
ceptibility. The n*^ order contribution to the matrix 11 
is by definition 

-(M(r,Q|fc)Xt(0,Q|g)/7(")) , (CI) 

where 

f/(n) = y_J_ rr / dT^H'^{T^) . (C2) 

Out of all possible contractions in (CI) we want to select 
only those that give rise to bubble diagrams. We denote 
their contribution by n|^^^. We wish to show that 

nfPA(r,g|^,g) = (KoKo..-oKono)(r,Q|^,g), (C3) 

where K is defined as 
Ka^{T, Q\k, q) = Yl n^^(r, Q\k, k')V^p{Q\k\ g), (C4) 

k' 

and o denotes a convolution and simultaneous matrix 
multiplication 



{KoK)^p (T,Q\k,q) = ^Y. j dnK„^{T - TuQ\k,k') K^p{TuQ\k' ,q). 



(C5) 



Fourier transforming (C3) then gives 

. X The basic identity underlying the inductive proof of 

Uf^{^{iLOn,Q\k,q) = (k*K*- ■ ■*K*U°YiUn,Q\k,q), (C3) is 



(C6) 



where the convolution * is defined as 

{K ^K)^^(iuJn,Q\k,q) = ^^Kcy^(iujn,Q\k,k') 



X{r,Q\k)U 



in) 



RPA 



{KoKo---oKoX) {r,Q\k), 

(Cll) 



Summing over n then leads to 



K^j3{iujn-,Q\k' ^q). (C7) where the contraction notation indicates that only con- 
tractions compatible with our RPA-like summation have 
been carried out. Eqn (C3) is clearly a direct consequence 
00 of (Cll). We now prove (Cll) by induction. For n = 1 

Il^^^{iuOn^Q\k^q) = Uf^^-^{iu;ri^ Ql^^q)- (C8) we prove by a lengthy but straightforward calculation 

that 



n=0 



This can be written as 



nRPA = ^ * = (/ _ K)-' * 



X{r,Q\k)U^ 



1) 



where 



(C9) 



(CIO) 



{KoX){T,Q\k). (C12) 



RPA 



The induction step is then straightforward. We have 



X(r,Q|/c)/7(^+^) 



RPA 



J dTn 



+ 1^^4(^/1+1) 



RPA 

(C13) 
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where 



F = X(r,Q|fc)C/(") 



(C14) 



RPA 



The combinatorial factor n + 1 cancels exactly against 
the l/(n + 1) in the definition of /7(^+^). Now using; the 
induction assumption (Cll) in (C14), the final contrac- 
tion in (C13) reduces to the induction start n = 1, thus 
establishing the validity of (Cll) for n + 1. 



1 - n{Ek) - n{Ek+Q) 



n{Ek) - n{Ek+Q) 

ILUn — Ek -\- Ek-\-Q ' 



(D2) 



(D3) 



APPENDIX D: ELEMENTS OF Us 

Introducing the definitions 

n{Ek) ^n{Ek+Q) - 1 



n{Ek^Q) - n{Ek) 

iuJn -\- Ek — Ek-\-Q 



B__{k,k^Q) 



Eh — E-^ 



k+Q 



(Dl) the explicit elements of the matrix, Eq. (89), are 
I 



Iif^{iuJn.Q\k,k') 

Iif^{iun,Q\k,k') 
IifS^n.Q\k,k') 

nfi(icjn,Q|/^,/^0 = - 

Ii^^{iujn,Q\k,k') = 

nf3(zCJn,Q|^,^0 = 

nfi(icjn,Q|^,^0 

Ui^{iu;n,Q\k,k') 

Ul^{iUn,Q\k,k') 



^ Zk^k^Q Baa' {k, k^Q) {8k,k' - h-k'-o) , 

cr,cr' = ± 

cF\k>Zl>j^Q Baa'ik' ,k' ^ Q){5k^k' - h-k'-o), 



- ^ aXkZ^j^Q Baa'{k,k ^ Q)[Sk,k' -h-k'-q) , 

{Z^j^qZ'"' 8k^k' + crcr' XkXk+qh-k' -q) B^a' (k, k + Q), 
ct^A/c+qZ^"^ Baa' (k, k^Q) [Sk,k' - h-k'-o) , 



(D4) 



a^a' — ^ 



crcr'Afe+QA/c Bcra'{k,k ^ Q)[Sk,k' -h-k'-o) , 

=± 

■ ^ Cr'Afe'+Q^^/"^ (^^ + Q) {h,k' - h-k'-^ 

a,a'^± 

y^ ^A^'^^fe+Q k^Q) {Sk,k' - Sk-k'-q) . 

I 



a,a' — :h 



a,a' — :L 



(D5) 
(D6) 
(D7) 

(D8) 
(D9) 
(DIO) 

(Dll) 
(D12) 
(D13) 



APPENDIX E: FURTHER CONTRIBUTIONS TO 
THE XXZ SPIN CHAIN RESPONSE 



The calculation described in the main text leads to 
a response that features sharp thresholds even at finite 
temperature. Physically one expects these thresholds 
to be absent at T 7^ 0. The diagrams we consider 
are incapable of capturing this effect because they in- 
clude poles that only depend on two single particle en- 



ergies. This limits the response at positive frequencies 
to the regions m.m{Ek Eq) < uo < max(£^/c + Eq) and 
< < max(£^A: — Eq)^ for general k^q. Coupling to 
decay channels involving more than two particles should 
alleviate this problem. 



One can evaluate the two-loop self energy correction 
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to the propagator. This includes diagrams of the form 
In analogy with the one-loop calculation, we define a self 



Here we have used the boson occupation factor, nB(e) = 
l/(exp(/3e) — 1). The remaining elements of the two loop 
self energy are obtained via 

= (E4) 
Sl^ = (Si2)*. (E5) 

These contributions have poles that feature three single 
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